rm(list=ls())

## install_version("devtools", version = "1.13.3", repos = "http://cran.us.r-project.org")
## require(devtools)
## install_version("msm", version = "1.6.4", repos = "http://cran.us.r-project.org")
## install_version("MASS", version = "7.3-47", repos = "http://cran.us.r-project.org")
## install_version("distr", version = "2.6.2", repos = "http://cran.us.r-project.org")
## install_version("truncnorm", version = "1.0-7", repos = "http://cran.us.r-project.org")
## install_version("VGAM", version = "1.0-4", repos = "http://cran.us.r-project.org")

library(sfaPA)
library(msm)
library(MASS)
library(distr)
library(truncnorm)
library(VGAM)
library(devtools)

##type is "both", "votes", or "words"
##CongressNum goes from 105-112
##shortrun does a short run, ~1-2 hours per model; set to FALSE to recover results in paper
##save.out saves the output or not

shortrun<-TRUE
Congress.use<-112
types.use<-"both"
save.out<-FALSE

## For a long run to reproduce the full results, uncomment the following
# types.use<- c("both","votes","words")
# Congress.use<-105:112
# save.out<-TRUE

for(type in types.use){
    for(CongressNum in Congress.use){
######################################################
######################################################
######### Load functions, set seed
######################################################
######################################################
        set.seed(1234)

        dubcent<-function(X){
            colmean.vec<-colMeans(X,na.rm=T)
            rowmean.vec<-rowMeans(X,na.rm=T)
            
            for(i in 1:dim(X)[1]) for (j in 1:dim(X)[2]) 
                X[i,j]<-X[i,j]-rowmean.vec[i]-colmean.vec[j]
            
            X-mean(X)
        }

######################################################
######################################################
######### Load in data
######################################################
######################################################

        load(paste("votemat",CongressNum,sep="_"))
        load(paste("dtm",CongressNum,sep="_"))

        if(type=="both"){
            dtm3 <- rbind(t(votes.mat),t(dtm2))
            missing.mat.use <- rbind(t(votes.mat==0.5),t(dtm2)*0)*1
        }

        if(type=="votes"){
            dtm3 <- rbind(t(votes.mat))
            missing.mat.use <- rbind(t(votes.mat==0.5))*1
        }

        if(type=="words"){
            dtm3 <- rbind(t(dtm2))
            missing.mat.use <- rbind(t(dtm2)*0)*1
        }

######################################################
######################################################
######### Format objects to pass to SFA
######################################################
######################################################


        unique.dtm<-sort(unique(as.vector(dtm2)))
        seq.vec<-c(-1:(max(dtm2)+2))
        for(i in 2:length(seq.vec)) ifelse(seq.vec[i]%in%unique.dtm, seq.vec[i]<-seq.vec[i-1]+1,seq.vec[i]<-seq.vec[i-1] )

        rank.mat<-0*dtm2
        for(i in 1:nrow(dtm2)) for(j in 1:ncol(dtm2)) rank.mat[i,j]<-seq.vec[dtm2[i,j]+2]
        empir.all<-0


######################################################
######################################################
######### Run SFA
######################################################
######################################################

        ##burnin: no. of burnin draws
        ##max.loop: no. of draws after burnin
        ##thin.save: thinning
        ##max.optim: use optim instead of HMC on part of model
        ##the rest: for internal use, will be updated for CRAN version

        if(!shortrun){
            set.seed(1234);obj.curr<- sparsemix(dtm3,missing.mat.use, max.loop=5000,
                                                burnin=3000, empir.cdf=empir.all,
                                                cutoff.seq=seq.vec, save.each=FALSE,
                                                thin.save=20, max.optim=500, seed=1234)
        }

        if(shortrun){
            
            set.seed(1234);obj.curr<- sparsemix(dtm3,missing.mat.use, max.loop=500,
                                                burnin=300, empir.cdf=empir.all,
                                                cutoff.seq=seq.vec, save.each=FALSE,
                                                thin.save=10, max.optim=100, seed=1234)
        }

######################################################
######################################################
######### Save object
######################################################
######################################################

        file.save<-paste("obj",type,CongressNum,sep="_")
       if(save.out) save(obj.curr,file=file.save)


    }} #Close out word type and senate loops

session_info()
